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SUMMARY 

A finite-element alternating method is presented for two-dimensional 
mode -I crack problems. An analytical solution for an arbitrary polynomial 
normal pressure distribution applied to the crack faces is obtained and used 
as the basic solution in the method. The method is applied to several crack 
problems to study its efficiency and the results are compared to accurate 
stress -intensity factor solutions in the literature. The method gave 
reasonably accurate stress -intensity factors and crack opening 
displacements with minimal computing effort. Because the method must model 
only the uncracked body, finite-element models with many degrees of freedom 
are not warranted and therefore, the method has been implemented on 
personal computers. 


INTRODUCTION 


Damage tolerant design concepts for aerospace structural components 


are widely used to avoid catastrophic failures during service. Stress- 
intensity factors are important parameters in these designs, because they 
are necessary to predict fatigue crack growth rates and fracture strengths. 

In two-dimensional analyses, several methods are available to 
calculate the stress -intensity factors of cracked components. Several 
stress- intensity factor compendia [1-4] are also available. However, there 
are always situations where the stress -intensity factor for a particular 
crack configuration and loading is not readily available and some 
approximations and estimates need to be made. Therefore, the search for 
algorithms, new methods and computer programs that are fast, accurate and 
efficient continues. 

Recent literature shows that the finite-element alternating method 
(FEAM) is a powerful, accurate and computationally efficient [5-8] method 
for three-dimensional (3-D) analysis. This method was successfully applied 
to embedded-, surface-, and corner -cracked solids with elliptic or part- 
elliptic cracks. The method is based on the Schwartz -Neumann alternating 
method [9,10], a numerical method that alternates between two solutions to 
satisfy the required boundary conditions. Solution 1 is usually a 
continuum solution for a cracked infinite plate or solid. Solution 2 uses a 
numerical method such as the finite element method to model and analyze the 
uncracked plate or solid with the same configuration. The method alternates 
between these solutions to satisfy the required boundary conditions of the 
original problem. The 3-D FEAM was shown to be very efficient from the 
modeling point of view and also gives accurate stress- intensity factors [4- 
8]. The method is also more economical than the conventional finite-element 
method. Judging from the performance of the 3-D FEAM, a two-dimensional (2- 
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D) version of the method appears to be desirable. The purpose of the paper 
is to develop a 2-D version of FEAM, test the feasibility of using the 
method on personal computers, and study its efficiency by applying it to 
various crack configurations under mode-1 loading. 

First, the analytical solution for a crack in an infinite plate 
subjected to arbitrary crack face pressures is obtained. Next, the details 
of the FEAM are presented. Third, FEAM is applied to two-dimensional crack 
configurations for which exact or accurate solutions are available, to study 
the accuracy of the stress -intensity factors produced by the method. 

Several computational aspects of the method are also discussed. 
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LIST OF SYMBOLS 
half-length or length of crack 
normalized stress- intensity factor 
half-height of plate 

integral of residual norm defined in Equation (25) 
mode -I stress- intensity factor 

length of j c ^ element along region of crack or external 
boundary 

pressure distribution applied to crack faces 

two symmetric concentrated loads on crack faces at a 

distance s from center of crack 

number of elements along crack boundary 

direction cosines of normal to boundaries along x- and 

y- directions, respectively 

remote uniform applied stress 
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u,v 

W 

x.y 

z 

Matrices 

(A) 

[K] 

{Q^} 

(uJ) 

Greek Letters 
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*1 


[O) 

v 

Superscripts 

T 

j 


displacements in x- and y- directions, respectively 
half -width or width of plate 
Cartesian coordinates 

complex variable, z- x + iy, i -(-l)^/ 2 

coefficients of polynomial residual pressure distributions 
assembled stiffness matrix of finite-element model 
assembled load vector for iteration 
finite-element displacement vector for iteration 

crack opening displacement at center of 
crack, (z-0) . 

maximum normalized crack opening displacement, 

Aj - E*/[4(l-t/ 2 )Sa] 

Cartesian stresses { a x Oy a X y } T 

Westergaard stress function, “ ft dz . yb 9 ■ dit 

dz 

shear modulus 
Poisson's ratio 

denotes transpose of a matrix 
jt* 1 iteration 


ANALYTICAL SOLUTION 

Consider an infinite plate with a crack of length 2a as shown in 
Figure 1. Use of the alternating method requires the solution for stresses 
and displacements everywhere in the plate due to an arbitrary polynomial 
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pressure distribution, p y , applied to the crack faces. The applied 
polynomial pressure distribution is assumed to be of the form 
N 

p y - S A n (x/a) n ( 1 ) 

n-0 

To obtain the stresses and displacements everywhere in the infinite 
plate due to this arbitrary pressure distribution, the analytical solution 
for a typical term, (x/a) n , must be obtained. To do this, consider the 
solution for a pair of normal forces P acting on the crack faces at a 
distance s from the center of the crack as shown in Figure 2. The 
Westergaard stress function, \tp,is [1] 




P < a 2-.2)l/2 

(z-s) (z 2 -a 2 ) A /2 


( 2 ) 


where z - x + iy. 

Using Vp as a Green's function, the stress function 
arbitrary pressure distribution of Equation (1) is given by 

a 




ir a* 


1 

— 2 — 7 


-jj2 J * s0 ( a2 ' s2 ) 1//2 )/(z-s) ds 


-a 

The singular integral in Equation (3) is of the form 




for an 


(3) 




a 

(s n (a 2 - s 2 ) 1/2 )/(s-z) ds 
-a 



(4) 


Integration over a contour enclosing the crack gives 
I n (z) - ir [ z n (z 2 -a 2 )V2 . Gm(z) ] 
where G^iz) is the principal part of { z n+ ^ (1 -a 2 /z 2 ) V 2 


(5) 

at z - ® [ 11] . 
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The Gco(z) in Equation (5) can be conveniently expressed for odd and 
even values of n as 


m+1 


Ga>(z) - z 2m+2 2 C k ( -a 2 /z 2 ) k / k! 

k-0 

when n - 2m+l, 


( 6 ) 


m 


G«o(z) - z 2m+1 2 C k (-a 2 /z 2 ) k / k! 

k-0 

when n - 2m, and 


(7) 


{ 1 for k-0, 

( 8 ) 

(1/2) (1/2-1) (1/2-2) ... (1/2-k+l) fork-1,2,3... 


The stress function in Equation (3) for any integer power n can 
then be written as 

- (z/a) n + G c0 °(z)/[ a n (z 2 - a 2 ) X /2 ], for odd n 


tf(z) 


1 


(z/a) n + G 00 e (z)/[ a 11 (z 2 - a 2 )^/ 2 ], for even n. 


(9) 


The stresses at any location z in the plate are given by 
a x - Re (i>) - y Im(^' ) 

Oy - Re (\f>) + y Im(ifr') (10) 

<7 xy - - y Re(tf') 

where Re( ) and Im( ) denote the real and imaginary parts, respectively, of 
the function in the parentheses, and 

V ~ ( 11 ) 

dz 

The stress -intensity factors are given by 

K - (2*) 1 / 2 Lim (z-a) 1 / 2 ^(z) (12a) 

z*~a 

for the crack tip at x - a, and 
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(12b) 


K - (2ir ) X / 2 Lira (-i)( z+a) 1 / 2 tf(z) 
z«^-a 

for the crack tip at x - -a. 

The stress- intensity factors for each of the distributions (x/a) n 
have been computed for n-1 to 8, and are presented in Table 1. These 
results agree identically with those obtained by Isida [11] and given in 
reference 4 (page 189) for n~ 0 to 3. 

The displacements at any point z can also be obtained from the 
Westergaard stress function as 

2/iu - (l-2i/) Re (VO - y Im (VO 

2/iV - 2(l-i/) Im (VO - y Re (VO (13) 

for plane strain. In Equation (13) is the shear modulus, u is the 
Poisson's ratio and 
V> — J $ dz 

For plane stress conditions, i/ is replaced by i//(l+i/) in Equation (13). 

The integrations involved in Equation (13) are straightforward and 
explicit expressions are not presented here. The maximum crack opening 
displacements, 5, at z«0 were obtained for each of the terms (x/a) n and are 
presented in Table 1. 


FINITE -ELEMENT ALTERNATING METHOD 

As previously mentioned, this method alternates between two solutions to 
satisfy the boundary conditions of the problem. Solution 1 is for an 
infinite plate with a crack subjected to arbitrary pressure distribution 
(Equations 9 and 10). Here for solution 2 the finite-element method is 
used. The procedure that is followed by the FEAM is given in reference 5 
and is summarized here and in Figure 2 for completeness. 
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Step 1 . Analyze the same configuration and loading as in the given 
problem, but without the crack , using the finite -element method. 

Step 2 . The finite -element method gives the stresses everywhere in the 
plate including the line segment coincident with the crack in the original 
problem. Because only mode-I configurations are considered here, no shear 
stresses o X y are present on the crack line. 

Step 3 . If the tractions on the crack faces are negligibly small (i.e. 
smaller than a prescribed tolerance level) , stop and calculate the sum of 
the stress -intensity factors computed so far. If the tractions are not 
negligible go to step 4. 

Step 4 . To free the tractions on the crack faces, the crack face normal 
tractions, py, computed in step 3, must be erased. To do this, the negative 
of crack- face normal tractions (that is, p y R - -p y ) are applied in the 
analytical solution. The tractions p y R are expressed in polynomial form as 
N 

Py R (x) - Z A n (x/a) n - (P) T (A) 
n-0 

where (P) T - {1 (x/a) (x/a) 2 . . . (x/a) N } T (14) 


The coefficients (A) in Equation (14) are then calculated using a 
least squares procedure by 

(A) - [H]' 1 (D) (15) 

where 

N c r 

[H] - Z /{ P) T {P) dx (16) 

j-1 J 


ID) 


N c r 

- Z / (P) T Py R (x) dX 

j-1 J 
'j 


(17) 
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where N c - number of elements and ij is the length of the element along 

the region of the crack. The integrals involved in Equations 16 and 17 are 
evaluated by Gaussian quadrature because discrete values of py R (x) are 
available in each of the N c elements. 

Step 5 . The coefficients (A) determined in Equation (15) are then used to 
calculate the stress-intensity factors for the current (i^) iteration, K^, 
from 

N 

K 1 - 2 (k w ) n A n l (18) 

n-1 

where (k w ) n are the stress - intensity factor weights given in Table 1 for 

each of the polynomial functions (x/a) n . 

Step 6 . The crack face normal tractions Py^ in step A create tractions on 

all the boundaries of the region of interest. These tractions are 

calculated for each of the polynomial functions by using Equations 9 and 10, 

at any point z on the boundary as 

(o) - [M] (A) (19) 

(3,1) (3 ,N) (N, 1) 

where 

(o) T - {a x o y 0 xy ) T (20) 

and [M] are the tractions at any point z on the boundary due to unit values 
of each of the polynomial pressures. The normal (a n ) and tangential (o t ) 
tractions at any point on the boundaries of the region of interest are 
calculated by using 

{o n a t ) T - [q] (a) (21) 

where 
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( 22 ) 


n x 

0 

n y 1 

0 

n y 

n x 


where n x , riy are the direction cosines of the normal to the boundary with 
respect to the x and y axes, respectively. 

Step 7 . To satisfy the traction free boundary conditions the tractions 
created from step 6 by the residual pressures p<y^ on the crack face need to 
be erased. Therefore, the negative of these tractions are considered as the 
prescribed tractions for the finite-element model of the uncracked plate. 
They are conveniently expressed in terms of the nodal forces on each of the 
elements on the boundaries of the finite-element model by 

{ Q ) j - ~[G]j (A) (23) 

where 


[GJi - 


/ 


[N] T [q] [M] ds 


(24) 


where the subscript j refers to the j 1 -* 1 element, [N] are the element shape 
functions, [q] is the direction cosine matrix in Equation (22), [M] is 
defined in Equation (19), and ij is the length of the side of element j on 
the boundary. 

The nodal forces {Q)j on the j element on the external boundaries 
are treated as applied forces and are assembled to form a new load vector 
for the finite-element model of the uncracked body. Using this load vector, 
the uncracked body is analyzed once again (step 1). This is the start of 
the next iteration. The iterative process is continued until the crack face 
tractions in step 3 are negligibly small. In the converged solution, the 
stress-intensity factors are simply the sum of the stress-intensity factors 
for all iterations. 

The criteria for convergence of the FEAM were formulated in terms of 


I 
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the total magnitude of the crack face pressures and the incremental 
contributions to the stress -intensity factors. Define the integral of the 


residual normal tractions py R in the region of the crack i c by 


IR 1 


/' ^ 11 


dx 


(25) 


^c 

for the i th iteration. If the integral lR i+ ^ is greater than Ir* and Ir*" 1 
the algorithm is terminated because of nonconvergence. The algorithm is 
continued as long as 


lR i+1 < V > IR 1 ’ 1 
lR i+1 < IR 1 < IR 1 ' 1 

and is terminated when (Ir/Ir^) is less than a predetermined tolerance. 

Because the stress- intensity factors are sought by this method, the 
ratio of the absolute value of the stress- intensity factor increment from 
the current iteration i, to the sum of the stress - intensity factors computed 
so far can be used in a criterion for convergence as 


i- 1 

r - K* / 2 KJ <0.01, for example. (26) 

j-1 

Numerical experimentation showed that ( Ir/Ir^" ) is a better 
criterion than Equation (26) for monitoring the convergence and terminating 
the algorithm. This is because a choice of r-0.01 does not mean that the 
final stress-intensity factors are within one percent of the "true" stress- 
intensity factors. Therefore, the ratio (Ir/Ir*) is used and when the ratio 
(Ir/Ir 1) - 0.01, the algorithm is terminated. The dependence of the final 
stress -intensity factors on the value of the ratio (Ir/Ir*) used to 
terminate the algorithm is briefly examined later in the paper. 

Obviously the stress- intensity factor solutions obtained by the 
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FEAM depends on the finite element mesh refinement for the uncracked plate. 
To ensure that the FEAM converges to the correct solution, at least two 
finite element models, coarse and fine, of the uncracked plate are 
desirable. This aspect also is examined in a later section of the paper. 

Computational Aspects of the FEAM 

In the finite element part of the method, the 8-noded isoparametric 
parabolic elements were used to model the uncracked body. Some 
computational aspects that make the method efficient are briefly discussed. 

Decomposition of the Global Stiffness Matrix .- In the FEAM the 
uncracked body is analyzed by the finite -element method several times (the 
number of iterations required for convergence of stress - intensity factors) 
with a different load vector each time. Therefore, the finite -element 
stiffness equations can be written as 

[K] {u°, u 1 , u 2 , . . . ) - {Q°, Q 1 , Q 2 , . . . ) (27) 

where [K] is the assembled stiffness matrix, and {u*-} is the displacement 
vector corresponding to the load vector (Q^) for the i*^ 1 iteration. 

Equation (27) can be very efficiently solved by finding the Cholesky 
factors of the stiffness matrix [K] soon after the stiffness matrix is 
assembled. This decomposition needs to be performed only once. Thereafter, 
the displacements { u M for the load vector (Q*) can be obtained by back 
substitutions. The displacements (u*} can be obtained very inexpensively 
because the back substitutions consume far less CPU Time than the 
decomposition [4-8] . 

Boundary Tractions .- In step 7 of the method the nodal forces on each 
element on the boundaries are obtained. This computation can be 
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efficiently performed by calculating the [G] matrices in Equation (24) for 
all elements on the boundaries before the start of the first iteration. In 
each iteration only one matrix multiplication in Equation (23) is required 
for each element on the boundaries. 

Stress Smoothing . - In step 2 of the method the distribution of py along the 
crack line is calculated from the finite-element solution. The nodal 
stresses computed in a finite -element solution can be inaccurate. 

Therefore, the stresses are computed at the 2x2 Gaussian points in an 
element and are extrapolated to the element nodes. ( See references 12 and 
13 for details on stress smoothing.) 

Fictitious Pressures . - In the case of edge cracks, only one crack tip 
exists. The analytical solution used in the FEAM is based on an interior 
crack in an infinite plate. Therefore, to use the method for edge cracks, 
the crack face pressures were fitted onto the fictitious part of the crack 
which lies outside the analysis region. A simple linear extrapolation with 
Py^ values at x - 0 and x - a was used to obtain the values for -a < x < 0, 
as shown in Figure 3(a). This method is preferred over the extrapolation 
from a least-squares fit because extrapolation outside the analysis region 
is known to be inaccurate [5]. 

This concept of the fit in the region -a < x < 0 was also used from 
the view of computational convenience for crack configurations that are 
symmetric about the y - 0 line. For these configurations, only the even 
powers of n in the polynomial fit of Equation (14) for py^ need to be 
retained. However, it is convenient from a programming standpoint to retain 
the complete degree polynomial for Py^, and to use the fictitious 

pressures in the region -a < x < 0 that are symmetric to those in the 
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region 0 < x < a, as shown in Figure 3(b). By using this type of a fit, the 
coefficients A n corresponding to the odd powers of x/a in Equation (14) 
automatically vanish, leaving only the even powers. 

Crack Tin Location .- Because the finite element method analyzes the 
uncracked plate, the location of the crack tip needs to addressed. The 
simplest way to define the crack tip is to assume that the crack- tip 
coincides with an end node in the finite element model as shown in Figure 
4(a). In this case the shaded elements are used to evaluate the tractions 
Py R in the crack region for a crack length of a^ ( see Figure 4(a) ). The 
requirement that the crack tip coincide with a finite element node is, 
however, unnecessary. Figure 4(b) shows a slightly longer crack than in 
Figure 4(a) t one whose tip does not coincide with a finite element node. In 
this situation one additional element needs to be used in the residual 
stress fit as shown by the shaded elements in Figure 4(b) . This technique 
is similar to that used in the 3-D case in reference 8 where the elements on 
the crack plane do not conform to the exact shape of the elliptical or 
part-elliptical crack front. 


EVALUATION OF THE FEAM 

To evaluate the effectiveness of the finite-element-alternating 
method, the method is applied to a variety of mode -I problems for which 
accurate solutions are available. First, the method is applied to three 
commonly encountered fracture specimens with rectilinear boundaries. In 
this application the sensitivity of the solution to finite-element mesh 
refinement and the degree of polynomial fit to the residual pressures is 
studied. Next, the method is applied to a C-shaped edge cracked specimen to 
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evaluate its effectiveness when curved boundaries are involved. Third, the 
method is applied to problems involving stress concentrations. Last, some 
observations on the computational efficiency of the method are made. 

In all cases plane strain conditions were assumed and a Poisson's 
ratio of 0.3 was used. 

Specimens with Rectilinear Boundaries .- Figure 5 shows three common 
fracture specimens with rectilinear boundaries: center -cracked tension 
(CCT), single-edge cracked tension (SECT), and three-point bend specimen 
(TPB). The FEAM was applied to these specimens to evaluate the method. 
Because of symmetries one quarter of the CCT specimen and one half of the 
SECT and TPB specimens were considered in the analysis. Figure 6 shows both 
coarse and fine finite-element idealizations. The coarse mesh had 65 nodes 
and 16 elements; the fine mesh had 173 nodes and 48 elements. These two 
meshes were used to analyze both the CCT and SECT specimens with a crack 
with a/W ratio of 0.5. 

Because of symmetry the v- displacements of all nodes along the y - 0 
line were prescribed to be zero for both CCT and SECT specimens. For the 
CCT specimen the u- displacements of all nodes along the x - 0 line were 
prescribed to be zero. For the SECT specimen only one node is prescribed 
to have zero u- displacement to prevent rigid body motions. The boundaries 
that must be made traction free are: (1) x - W and (2) y - H for the CCT 
specimen and (1) x -0, (2) x - W and (3) y - H for the SECT specimen. 
Therefore the [G] matrices of Equation (24) for the elements along these 
lines were calculated before the start of the iterative process. 

Table 2 presents the normalized stress- intensity factors, 
K/[S(jra)l/^] ( obtained with the two meshes and for different degrees of 
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polynomial representation, N, in Equation (14). For the CCT specimen 
excellent accuracy was obtained even with only a quadratic representation 
(N-2) and the coarse (16 -element) mesh. Furthermore, convergence to this 
accuracy, (Ir/Ir^) < 0.01, was obtained with only 4 iterations. Figure 7 
shows the convergence of the residual pressure for both CCT and SECT 
specimens. For the SECT specimen, on the other hand, a very good solution 
was obtained with N - 4 for the coarse mesh and with N - 5 for the fine 
mesh. The convergence to this accuracy was achieved with 11 iterations, 
which was much slower than with the CCT specimen. This was to be expected 
because the SECT specimen has an additional traction- free boundary (x - 0 
line). Because accurate solutions for both specimens were produced with a 
polynomial of the fifth degree, N— 5 was used in the rest of the paper. 

The use of a single fine mesh has a distinct advantage with this 
method. After the global stiffness matrix of the uncracked body is 
assembled and decomposed, a series of crack lengths can be analyzed in a 
single run. This is economical because the decomposition of the global 
stiffness matrix is the most expensive part of this method [4-8], For each 
crack length, however, the [G] matrices of Equation (24) need to be 
computed. This procedure was followed to obtain stress -intensity factors 
for 6 crack lengths in a single run, for each of the CCT, SECT, and TPB 
specimens. The normalized stress -intensity factors and maximum crack- 
opening displacement (COD) for various a/W ratios are presented in Tables 
3 and 4 for the CCT and SECT specimens, respectively. Tables 3 and 6 present 
the corresponding results for the TPB specimen. Reference results from the 
literature are also included in these tables for comparison. Excellent 
agreement was obtained for each of the five a/W ratios for each of the 
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configurations except for long cracks with an a/W ratio of 0.75. The 
maximum percent difference is observed to be about 12.5 percent for the 
stress -intensity factor. When the crack is long, a much finer mesh near 
the crack tip than that of Figure 6(b) is probably needed to obtain an 
accurate solution. 

To ascertain that this is the case a new finite element model of the 
uncracked plate was developed and is shown in Figure 6(c). This model had 
173 nodes and 48 elements as in Figure 6(b) but the mesh spacing is non- 
uniform in the x-direction. The stress-intensity factors were recomputed 
for the SECT specimen with this model for a/W ratios of 0.5, 0.625, and 
0.75. Table 3 presents these results parenthetically. The new value of the 
normalized stress- intensity factor for a/W - 0.75 is now in excellent 
agreement with the reference value from the literature. The values for 
a/W-0.5 and 0.625 are about 3 percent lower than the reference value. This 
is expected because in the non-uniform model for a/W ratios of 0.5 and 
0.625 the residual pressures were fit with stress values from 2 and 3 
elements, respectively, while in the uniform model the fit was made from 4 
and 5 elements. Thus, for a/W ratios up to 0.6, the uniform model performed 
well and for a/W ratio of 0.75 the non-uniform model gave accurate results. 
Therefore, a judicious choice of the uncracked finite element models is 
required to produce accurate stress-intensity factors with FEAM. 

Recall that the results in Table 3 were obtained with an (Ir/Ir^) 
ratio of 0.01. To study the effect of the termination value of (Ir/Ir 1 ) , a 
more stringent value of 0.001 was used. The stress- intensity factors for 
the SECT specimen with a/W ratios of 0.5, 0.625, and 0.75 ( with the uniform 
mesh , Figure 6(b) ), were calculated to be 2.845, 4.620, and 9.536, 
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respectively. These values differed by less than one percent from the 


values given in Table 3. Therefore, a value of (Ir/1r^) ^ 0.01 was 
sufficient to produce converged values of stress -intensity factors and is 
used to terminate the algorithm in the rest of the paper. 

Specimen with Curved Boundaries .- Figure 8(a) shows a C-shaped edge cracked 

specimen proposed by Kapp et. al [15]. This specimen was analyzed to 

evaluate the applicability of the method to problems involving curved 

boundaries. Figure 8(b) shows a finite element model of the uncracked body. 

Because of symmetries only the top half of the specimen was modeled. 

Because the pin holes and the region to the left of the pin are small and 

remote from the crack tip, they were not modeled. Figure 8(b) shows the 

finite element model which had 193 nodes and 54 elements. The curved 

boundaries in the specimen were modeled by piecewise linear segments. 

The v- displacements of all nodes on the y - 0 line were prescribed 

zero. One node was constrained in the u-direction to prevent rigid body 

motion. A single concentrated load P was applied at node E of the model 

(see Figure 8(b)). The [G] matrices of Equation (24) were calculated for 

all elements along the curved boundaries BC and AD and along the straight 

boundary CD. As before the global stiffness matrix was assembled and 

stress - intensity factors for 6 normalized crack lengths ranging from a/W - 

0.125 to 0.625 were calculated in a single run. Table 7 presents the 

normalized stress-intensity factors, K/K(a) , where [15] 

K(a) - P { 3(X/W) + 1.9 + 1.1a } [1 + 0.25 (1-a) 2 (l- ri /r 2 )] 

W 1 / 2 

with 

a - a/W ( See Figure 8(a)). 

The present results are generally lower than the results of Kapp et. 
al [15] by about 4 percent. The accuracy of the results in ref. 15 is given 
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as ±1.5 percent. In the present analysis the pin holes and a small portion 
of the specimen were not modeled and the curved boundaries were modeled as 
piecewise linear, which may account for these small differences. The 
present results, however, indicate the applicability of the method to 
specimens with curved boundaries. 

Specimens with Stress Concentrations .- To evaluate the FEAM in the presence 
of the additional complication of stress concentrations, two additional 
configurations have been analyzed ( see Figure 9) . The first configuration 
is that of two symmetric cracks emanating from a circular hole in a finite 
width plate subjected to remote tensile loading (Figure 9(a)). The second 
configuration is that of a crack emanating from a semi-circular notch in a 
finite width plate ( Figure 9(b)) subjected to remote tensile loading. 

Figure 10 presents the finite -element idealization of one quarter of 
the plate for the symmetric crack configuration and one -half of the plate 
for the semi-circular notch configuration. The finite-element model had 213 
nodes and 60 elements. For both configurations the v-displacements on the 
y-0 line (AB in Figure 10) Were prescribed to be zero. For the cracks from 
a hole, the u-displacements at all nodes on the line x - -R (line ED in 
Figure 10) were prescribed to be zero. For the crack from a notch, only one 
node was prescribed to have a zero u-displacement to prevent rigid body 
motions . 

In both configurations the lines BC, CD, and EA had to be made 
stress free. Hence, the [G] matrices of all elements along these lines were 
calculated at the start of the iteration process. For the cracks from a 
circular hole, along the line DE symmetric boundary conditions require that 
the shear stresses be zero. To achieve this the [G] matrices of Equation 
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(24) for all elements along the line DE were calculated and the forces in 
the y- direction were used as the residual forces in Step 7 of the method. 

For the crack from a semi-circular notch, on the other hand, both x- and y- 
forces at all nodes along this line were used as the residual forces. 

The stress-intensity factors were calculated for several crack 
lengths for each of these configurations in a single computer run. Tables 8 
and 9 present the normalized stress-intensity factors for these 
configurations. They are compared with reference solutions from the 
literature [16,17]. Good agreement was obtained between the two sets of 
results. These results suggest that the FEAM can give accurate stress- 
intensity factors for problems involving stress-concentrations. 

In all the examples analyzed so far the crack- tip in the original 
problem was always assumed to end at an end-node in the finite element model 
(see Figure 4(a) ). As previously mentioned, this restriction can be 
removed by the technique suggested in Figure 4(b). To validate this 
technique, the problem of a crack from a semi-circular notch is analyzed 
with crack lengths that do not coincide with the finite-element nodal 
points. Table 10 presents the normalized stress- intensity factors obtained 
by this technique and compares them with those calculated by the boundary- 
force method [17]. When only one element is used to obtain py R tractions 
and when the crack region is contained within this one element, the stress- 
intensity factors are not very accurate. However, when more than one 
element is used to model the crack region the results are within one 
percent of those from reference 17. These results suggest that for any 
crack length one could obtain the stress - intensity factors using a single 
fine mesh idealization provided that more than one element is used to model 
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the crack region. 


Computer Program 

The FEAM algorithm described here was written using FORTRAN 77 and 
was developed on a Personal Computer (PC) using a FORTRAN Compiler. The 
program was executed on personal computers with 8086, 80286, and 80386 
microprocessors. Then the program was ported to a main frame computer with 
only one change in the source code. For the problem of two symmetric cracks 
from a circular hole ( Figure 9(a) ), the execution times on each of these 
computers for an analysis with six crack lengths ( finite-element mesh with 
213 nodes and 60 elements) are as follows. 

Computer PC- 8086 PC- 80286 PC- 80386 CYBER 850 

CPU Time ( Minutes ) 120 16 8 1.63 

These computing times demonstrate the versatility of the method and show 
that one could obtain accurate stress- intensity solutions on personal 
computers with this method. 

CONCLUDING REMARKS 

A finite -element alternating method (FEAM) for two-dimensional crack 
configurations is presented. The method is based on an analytical solution 
for infinite plate containing a crack subjected to arbitrary normal 
tractions. The complete analytical solution for the infinite plate is 
given. The details of the FEAM are explained along with some computational 
aspects that improve the performance of the method. 

The method is applied to several crack configurations for which 
reference stress -intensity factor solutions are available in the literature. 
Even though no attempt was made to optimize the finite element mesh, the 
method gave reasonably accurate stress -intensity factors and crack opening 
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displacements for all problems analyzed. A fifth- degree polynomial fit to 
the residual pressures was found to be sufficient. 

Stress- intensity factors for several crack lengths were obtained in 
a single computer run with one finite -element idealization. The method can 
give a nearly continuous distribution of stress -intensity factors as a 
function of the crack length in a single computer run. The method was 
developed and implemented on a microcomputer (personal computer) . The 
results demonstrate the feasibility of obtaining stress -intensity factors 
for a variety of cracked bodies on personal computers. 
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Table. 1: Stress -intensity factors, K, and the COD at the center of a 
crack in an infinite plate subjected to arbitrary polynomial 
pressure distributions, p y - (x/a) n . 


n 

K 

<^F72 

uS * 

2(l-i/)a 

0 

1 

i 

1 

± 1/2 

0 

2 

1/2 

1/6 

3 

± 3/8 

0 

4 

3/8 

3/40 

5 

± 5/16 

0 

6 

5/16 

5/112 

7 

± 35/128 

0 

8 

35/128 

35/1152 


The positive and negative signs in this table 
refer to the crack tips at x~ ± a respectively. 

The negative values are meaningful only in the presence of 
additional forces which prevent crack closure. 

Note that . ut ’ “ K/ f wal for odd values of n 

2(l-i/)a n+1 



Table 2: Normalized stress -Intensity factors, Fj , for the CCT 

and SECT specimens with coarse and fine meshes and with 
various degree of polynomial fits ( a/W-0.5; H/W-2.5 ). 

Fj - K/[S(wa )V 2 ] 


Degree of 
Polynomial 
fit, N 

Coarse 

Mesh 

Fine 

Mesh 

Reference 
Value [1-4] 

CCT Specimen 

2 

1 . 188 (4) a 

1.187(4) 


4 

1.188(4) 

1.187(4) 

1.189 

6 

1.188(4) 

1.187(4) 



SECT Specimen 


2 

2.964(13) 

3.027(13) 


3 

2.807(11) 

2.918(11) 


4 

2.808(11) 

2.994(12) 

2.828 

5 

2.752(11) 

2.831(11) 


6 

2.744(11) 

2.800(11) 



a Values in the parentheses are the number of iterations required for 
convergence to one percent accuracy in the stress-intensity factors. 



Table 3: Comparison of normalized stress -intensity factors, Fj, for the 
CCT and SECT specimens by the FEAM with accurate values from 
literature ( H/W-2.5 ). 

F i -K/[S(ira)V2] 


a/W 



0.125 

0.25 

0.375 

0.5 

0.625 

0.75 

CCT Specimen 

FEAM 

1.009 

1.039 

1.094 

1.187 

1.342 

1.626 

Reference 
Value [1-4] 

1.010 

1.040 

1.097 

1.189 

1.342 

1.617 

Percent 

Difference 3 

-0.01 

-0.10 

-0.27 

-0.17 

0.0 

0.56 


SECT Specimen 


FEAM 

1.223 

1.490 

1.971 

2.831 

4.587 

9.447 

Reference 




(2 . 75 7 ) b 

(4.349) 

(8.455) 

Value [1-4] 

1.221 

1.494 

1.975 

2.828 

4.481 

8.481 

Percent 

difference 

0.0 

-0.27 

-0.20 

0.11 

2.4 

11.4 


(-2.5) (-3.0) (-0.31) 


a The percent difference here is based on (FEAM Value - Reference 
Value)/ Reference value. 

k Values in the parenthesis were obtained with non-uniforra mesh. 
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Table 4: Comparison of maximum normalized crack opening displacement, 
Aj , for the CCT and SECT specimens by the FEAM with accurate 
values from literature ( H/W-2.5 ). 

A X - E5/[4(l-»/ 2 )Sa] 


a/W 



0.125 

0.25 

0.375 

0.5 

0.625 

0.75 

CCT Specimen 

FEAM 

1.009 

1.039 

1.093 

1.181 

1.324 

1.569 

Reference 
Value [1,4] 

1.009 

1.039 

1.094 

1.182 

1.322 

1.572 

Percent 

Difference 3 

0.0 

0.0 

-0.09 

-0.08 

0.15 

-0.19 

SECT Specimen 

FEAM 

1.546 

1.988 

2.922 

4.925 

9.893 

27.34 





(4 . 676) b 

(9.148) 

(23.94) 

Reference 
Value [1-4] 

1.605 

2.044 

2.986 

4.922 

9.649 

23.04 

Percent 

difference 

-3.7 

-2.7 

-2.1 

0.06 

2.5 

18.7 





(-5.0) 

(-5.2) 

(3.2) 


a The percent difference here is based on (FEAM Value - Reference 
Value)/ Reference value. 

b Values in the parentheses were obtained with non-uniform mesh. 


28 





Table 5: Comparison of normalized stress -intensity factors, Fj , for the 
TPB specimen by the FEAM with accurate values from 
literature ( H/W-2.5 ). 

Fl - K/KeP/WXira) 1 ^] 


a/W 



0.125 

0.25 

0.375 

0.5 

0.625 

0.75 

FEAM 

0.970 

0.992 

1.123 

1.408 

2.030 

3.759 

Reference 






(3 . 356) b 

Value [1-4] 

0.977 

1.007 

1.137 

1.416 

1.981 

3.349 

Percent 

Difference 3 

-0.72 

-1.5 

-1.2 

-0.56 

-2.5 

12.2 


( 0 . 21 ) 


a The percent difference here is based on (FEAM Value - Reference 
Value)/ Reference value. 

k Values in the parentheses were obtained with non-uniform mesh. 
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Table 6: Comparison of maximum normalized crack opening displacement, 

Aj, for the TPB specimen by the FEAM with accurate values from 
literature ( H/W-2.5 ). 

A X - E 6 / 4(l-i/ 2 )Pa(6/W) ] 


a/W 



0.125 

0.25 

0.375 

0.5 

0.625 

0.75 

FEAM 

1.344 

1.548 

2.009 

2.955 

5.128 

12.17 

(10.74) 

Reference 

Value [4] 

1.373 

1.594 

2.041 

2.946 

5.100 

10.115 

Percent 

Difference 3 

-2.1 

-2.9 

-1.6 

0.31 

0.55 

20.3 

(6.2) 


a The percent difference here is based on (FEAM Value - Reference 
Value)/ Reference value. 

k Values in the parentheses were obtained with non-uniform mesh. 
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Table 7: Comparison of normalized stress -intensity factors, F(oc) , 

for the C- shaped specimen by the FEAM with accurate values 
from the literature ( X/W - 0.5; ^/rj - 2; « - a/W ). 

F(<x) - K/{ P/(W 1 /2)(3 x/w +1.9 +l.l<x) [l+0.25(l-cc)2(i. ri / r2 )) , 


a/W 



0.125 

0.25 

0.375 

0.4375 

0.5 

0.625 

FEAM 

1.331 

1.861 

2.536 

3.005 

3.625 

5.741 

Reference 

Value [4,15] 

1.316 

1.942 

2.650 

3.120 

3.733 

5.777 

Percent 

Difference 3 

1.1 

-4.2 

-4.3 

-3.7 

-2.9 

-0.62 


a The percent difference here is based on (FEAM Value - Reference 
Value)/ Reference value. 
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Table 8: Comparison of normalized stress -intensity factors, Fj, 

for two symmetric cracks from a circular hole by the FEAM with 
accurate values from literature ( R/W - 0.25; c - R+a ). 

F T - K/[S(*c)V2] 


c/W 


F I 


FEAM 

Reference Value [16] 

Percent 

Difference 3 

0.3 

— 

1.078 

— 

0.3125 

1.088 

(1 . 109) b 

-1.9 

0.375 

1.203 

(1.195) 

0.7 

0.4 

— 

1.216 

— 

0.4375 

1.254 

(1.240) 

1.1 

0.5 

1.311 

1.285 

2.0 

0.5417 

1.355 

(1.330) 

1.9 

0.583 

1.405 

(1.380) 

1.8 

0.6 

— 

1.397 

— 

0.667 

1.528 

(1.515) 

0.86 


a The percent difference here is based on (FEAM Value - Reference 
Value)/ Reference value. 

k Values in the parentheses are interpolated values. 





Table 9: Comparison of normalized stress- intensity factors, Fj , 
for a crack at a semi-circular notch by the FEAM with 
accurate values from literature ( R/W - 0.25; c - R+a ). 

Fl “ K/tS^c) 1 / 2 ] 




Fl 


c/R 

FEAM 

Reference Value [17] 

Percent 

Difference 3 

1.2 

— 

1.526 

— 

1.25 

1.572 

(1.605) b 

-2.1 

1.3 

— 

1.695 

— 

1.5 

1.959 

1.965 

-0.3 

1.75 

2.336 

2.336 

0.0 

2.0 

2.836 

2 . 818(2 . 827) c 

0.6 

2.167 

3.273 

3 . 255 c 

0.6 

2.333 

3.831 

3 . 783 c 

1.2 

2.667 

5.549 

5 . 409 c 

2.6 

3.00 

9.102 

8 . 481 c 

7.3 


a The percent difference here is based on (FEAM Value - Reference 
Value)/ Reference value. 


b Values in the parentheses are interpolated values. 

c Values correspond to an edge crack in a plate without the notch. For 
c/R £ 2, the effect of the notch is insignificant. 
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Table 10: Comparison of normalized stress -Intensity factors, Fj, 
for a crack at a semi-circular notch by the FEAM when 
the crack tip does not coincide with finite-element nodes 
( R/W - 0.25; c - R+a ). 

F T - K/[S(wc)1/2] 


c/R 


*1 


FEAM 

Reference Value [17] 

Percent 

Difference® 

1 . 125(l) b 

1.413 

1.350 

4.6 

1.150(1) 

1.481 

1.407 

5.2 

1.200(1) 

1.551 

1.526 

1.66 

1.375(2) 

1.776 

1.790 

-0.77 

1.400(2) 

1.807 

1.825 

-0.96 


a The percent difference here is based on (FEAM Value - Reference 
Value)/ Reference value. 

k The values in the parentheses are the number of elements in the crack 
region. 
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Fig. 1. Crack in an infinite plate subjected to two normal concentrated 


Plate with a crack 
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Fig. 2: Procedure used in the alternating method. 



Fictitious 



h— a — •+*- a — h 


(a) Residual pressures for edge cracks. 



(b) Residual pressures for crack configurations 
that are symmetric about x = 0 line. 


Fig. 3: Fictitious pressure extrapolations. 
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Fig. 4: Elements used in determining the residual pressures for 
different length cracks. 



Fig. 5: Specimens with rectilinear boundaries. (H/W— 2.5) 
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Fig. 6: Finite -element idealizations used for specimens with rectilinear 
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C-shaped specimen and a finite- element idealization 
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Fig- 9: Crack problems with stress concentrations 
( R/W— 0.25, H/W-2.0 ) 
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